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Abstract 

In  High  Dose  Rate  (HDR)  brachytherapy  the  conventional  dose  optimization  algorithms  consider  the  multiple  objectives  inform  of  an 
aggregate  function  which  combines  individual  objectives  into  a  single  utility  value.  A.?  a  result,  the  optimization  problem  becomes  single 
objective,  prior  to  optimization.  Up  to  300  parameters  must  be  optimized  satisfying  objectives  which  are  often  competing.  We  use 
multiobjective  dose  optimization  methods  where  the  objectives  are  expressed  in  terms  of  quantities  derived  from  dose-volume  histograms 
or  in  terms  of  statistical  parameters  of  dose  distributions  from  a  small  number  of  sampling  points.  For  the  last  approach  we  compare  the 
optimization  results  of  evolutionary  multiobjective  algorithms  with  deterministic  optimization  methods.  The  deterministic  algorithms  are 
very  efficient  and  produce  the  best  results,  but  they  also  have  the  certain  limitations.  The  performance  of  the  multiobjective  evolutionary 
algorithms  is  improved  if  a  small  part  of  the  population  is  initialized  by  deterministic  algorithms. 


1.  Introduction 

High  dose  rate  brachytherapy  is  a  treatment  method  for 
cancer  where  empty  catheters  are  inserted  within  the  tumor 
volume.  Once  the  correct  position  of  these  catheters  is 
verified,  a  single  192Ir  source  is  moved  inside  the  catheters  at 
discrete  positions  (dwell  positions)  using  a  computer 
controlled  machine.  The  problem  that  we  consider  is  the 
determination  of  the  n  dwell  times  (which  sometimes  are 
called  as  well  dwell  position  weights  or  simply  weights)  for 
which  the  source  is  at  rest  and  delivers  radiation  at  each  of 
the  n  dwell  positions,  resulting  in  a  three-dimensional  dose 
distribution  which  fulfills  the  defined  quality  criteria.  In 
modern  brachytherapy,  the  dose  distribution  has  to  be 
evaluated  with  respect  to  the  irradiated  normal  tissues  and 
the  Planning  Target  Volume  (PTV)  which  includes  besides 
the  Gross  Tumor  Volume  (GTV)  an  additional  margin 
accounting  for  position  inaccuracies,  patient  movements,  etc. 
Additionally,  for  all  critical  structures,  either  located  within 
the  PTV  or  in  its  immediate  vicinity  or  otherwise  within  the 
body  contour,  the  dose  should  be  smaller  than  a  critical  dose 
Dcrit.  In  practice  it  is  difficult,  if  not  impossible  to  meet  all 
these  objectives.  Usually,  the  above  mentioned  objectives 
are  mathematically  quantified  separately,  using  different 
objective  functions  and  then  added  together  in  various 
proportions  to  define  the  overall  treatment  objective  function 
[1,2]. 

The  number  of  source  positions  varies  from  20  to  300.  It  is 
therefore  a  high  dimensional  problem  with  competing 
objectives.  The  use  of  a  single  weighted  sum  leads  to 
information  loss  and  is  not  generally  to  be  recommended, 
especially  for  non  convex  problems  and  for  those  cases 
where  objectives  have  not  the  same  dimensions  and  in 
addition  maybe  competing.  An  understanding  of  which 
objectives  are  competing  or  non-competing  is  valuable 
information.  We  therefore  use  multiobjective  evolutionary 
algorithms  in  HDR  brachytherapy.  One  algorithm  is  based 
on  the  optimization  of  dose-volume  histograms  (DVH), 
which  describes  the  distribution  of  the  dose  within  an  object, 
or  from  these  derived  distributions.  These  distributions  are 
evaluated  for  the  PTV,  the  surrounding  tissue  and  organs  at 
risk  from  a  set  of  up  to  100000  sampling  points  [3],  The 
calculation  of  the  DVH  requires  a  considerable  amount  of 
time  and  for  implants  with  300  sources  the  optimization 
requires  a  few  hours.  Another  limitation  of  this  method  is 
that  a  comparison  with  deterministic  algorithms  is  not 
possible.  We  have  therefore  considered  the  optimization  of 


the  dose  distribution  using  as  objectives  the  variance  of  the 
dose  distribution  on  the  PTV  surface  and  within  the  PTV 
obtained  from  a  set  of  150Ck-4000  sampling  points.  These 
functions  are  convex  and  a  unique  global  minimum  exists. 

In  the  past  comparisons  of  the  effectiveness  of  evolutionary 
algorithms  have  been  made  with  either  other  evolutionary 
algorithms  [4]  or  with  manually  optimized  plans  [1,  2],  We 
have  compared  the  Pareto  fronts  obtained  by  multiobjective 
evolutionary  algorithms  with  the  Pareto  fronts  obtained  by  a 
weighted  sum  approach  using  deterministic  optimization 
methods  [5], 

We  use  here  only  objectives  where  gradient  based 
algorithms  are  superior.  However,  we  must  consider  also 
critical  structures  partly  inside  the  target  or  close  to  it  which 
have  to  be  protected  by  excessive  radiation.  Other  objectives 
are  the  optimum  position  and  the  minimum  number  of  sources. 
In  such  cases  the  gradient  based  algorithms  can  not  be  used. 

2.  Methods 

2.1  Calculation  of  the  Dose  Rate 

The  dose  rate  around  each  of  the  small  cylindrical  shaped 
sources  is  dominated  by  the  1/r2  term  with  modifications  due 
to  absorption  and  scattering  in  the  surrounding  material.  The 
dose  value  d(r)  at  r=(x,  y,  z)  is: 

d(x)  =  |>,.K(r-r,.)  (1) 

In  (1)  (  is  the  position  of  the  kh  source  and  bfc  the  total 
number  of  sources.  K(r-q)  is  the  dosimetric  kernel  describing 
the  dose  rate  per  unit  source  strength  at  r  from  a  source 
positioned  at  q.  The  dwell  position  weight  w;  =  SkE  is 
proportional  to  the  strength  Sk  of  the  of  the  single  stepping 
source,  where  (  is  the  dwell  time  of  the  kh  source  dwell 
position  [6],  Because  of  the  high  dose  gradients  a  dose 
specification  at  a  single  point  inside  the  PTV  is  not  possible  in 
interstitial  brachytherapy.  For  this  reason  we  use  as  a  reference 
dose  Drefthe  average  dose  value  at  the  PTV  surface. 

2.2  Dose-Volume  Histogram  Based  Optimization  Using 
the  Conformal  Index 

A  conformal  Index  (COIN)  was  proposed  as  a  measure  of 
implant  quality  and  dose  specification  in  brachytherapy  [7], 
This  index  takes  into  account  patient  anatomy,  both  of  the 
tumor  and  normal  tissues  and  organs. 

We  describe  the  dependence  of  the  conformal  index  COIN 
on  the  choice  of  the  reference  dose  value  as  the  COIN 
distribution,  see  Fig.  1(b).  Usually  the  dose  values  are 
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normalized  to  Dref  and  are  given  either  as  fractions  or 
percentages  of  Dref. 

The  cumulative  dose  volume  histograms  of  the  PTV  and 
the  body  for  a  rib  implant  is  shown  in  Fig.  1(a).  Due  the 
rapid  decrease  of  the  dose-volume  histogram  (DVH)  of  the 
body  a  large  number  of  sampling  points  is  necessary  in  order 
to  calculate  with  a  high  accuracy  the  DVH,  the  COIN 
distribution  and  the  COIN  integral  at  dose  values  close  to  the 
reference  dose  value  and  above.  The  COIN  distribution  from 
the  DVHs  of  Fig.  1(a)  and  the  COIN  integral  is  shown  in 
Fig.  1(b). 


PTV 

BODY 


Fig.l.  (a)  Dose-volume 
histograms  of  the  PTV 
and  the  body  as  a 
function  of  dose.  (b)The 
corresponding  COIN 
distribution.  The  shaded 
area  to  the  right  of 
D/Dref=1.5  is  the  COIN 
Integral.  The  objectives 
are  maximum  COIN 
value  at  D=D,ef  and 
minimum  COIN  integral 
for  the  avoidance  of  high 
dose  values  in  the  PTV 
and  the  surrounding 
tissue. 


2.3  Dose  Statistics  Based  Optimization 

The  DVH  based  optimization  method  requires  a  large 
number  of  sampling  points  for  the  computation  of  the 
histograms  and  the  COIN  distribution  and  therefore  is 
computational  expensive.  We  have  developed  a  stratified 
sampling  approach  where  the  sampling  points  are  non 
uniform  distributed  and  which  reduces  the  number  of 
required  sampling  points  by  a  factor  of  5h-10.  Even  then  for 
implants  with  200t-300  sources  the  optimization  time  can 
reach  k-2  hours.  A  comparison  of  the  performance  with 
deterministic  and  gradient  based  algorithms  is  not  practical 
or  not  even  possible. 

Therefore  we  consider  another  set  of  two  objectives:  For 
the  conformity  objective  we  use  the  variance  fs  of  the  dose 
distribution  of  sampling  points  uniformly  distributed  on  the 
PTV.  In  order  to  avoid  excessive  high  dose  values  inside  the 
PTV  we  require  a  small  as  possible  dose  distribution 
variance  fy  inside  the  PTV.  Due  to  the  source  characteristics 
these  two  objectives  are  competing.  We  use  normalized 
variances  for  the  two  objectives: 


f  =  - 


7I  (d-m) 


(4) 

m2N*~!' 

Where  m  is  the  average  dose  value  and  N  the  corresponding 
number  of  sampling  points. 


optimization  mns  where  ws  varied  from  0  to  1  in  steps  of  0.05 
to  determine  the  shape  of  the  trade-off  curve.  A  problem  in 
using  deterministic  optimization  methods  is  that  the  solution 
contains  a  large  number  of  dwell  weights  with  negative  values. 
This  is  a  non  physical  solution.  In  the  past  either  constrained 
optimization  methods  were  used  or  a  correction  was  applied 
by  setting  to  0  all  negative  weights  in  each  optimization  step. 
A  constrained  optimization  method  increases  the  number  of 
parameters  by  a  factor  of  two.  The  correction  method  for  the 
negative  weights  reduces  the  quality  of  the  optimization 
results.  We  use  a  simple  technique  by  replacing  the  decision 
variables,  the  weights  v^,  with  the  parameters  w\  =  wk1/2 . 
Using  this  mapping  technique  we  avoid  non  feasible  solutions. 
For  this  unconstrained  optimization  we  use  the  Polak-Ribiere 
variant  of  Fletcher-Reeves  algorithm  or  the  Broyden-Fletcher- 
Goldfarb-Shanno  quasi-Newton  based  algorithm  [5].  These 
require  the  first  derivative  of  the  objective  function  with 
respect  to  the  decision  variables  to  be  calculated.  The 
derivative  of  the  normalized  variance  /  used  by  the  gradient 
based  optimization  methods  is: 


3  f 

3i vk  Nm 3  \ 


£(d,-m) 


3  d,  ,  dm 

m — -~d — - 
3i/vt  3 


(6) 


As  a  gradient  free  method  we  used  the  modified  Powell 
method  of  Numerical  Recipes  [5]. 


2.5  Multiobjective  Optimization  with  Evolutionary 
Algorithms 

The  population  of  multiobjective  algorithms  consists  of  strings 
storing  a  set  of  weights  for  each  source  dwell  position.  The 
weights  are  initially  produced  randomly  distributed  within  the 
interval  [0,1],  A  part  of  population  can  be  initialized  by 
solutions  of  the  deterministic  algorithms,  and  more  on  this  will 
be  written  further  [8]. 

In  our  algorithm  analysis  we  used  these  selection  mechanisms: 
-The  niched  Pareto  algorithm  (NPGA)  proposed  by  Horn  and 
Nafpliotis  [9]; 

-Strength  Evolutionary  Approach  algorithm  (SPEA)  by  Zitzler 
and  Thiele  [10], 

-Non-dominated  Ranking  Algorithm  by  Fonseca  and  Flaming 
(FFGA)  [11, 12] 

-Non-elitist  Non-dominated  Sorting  Genetic  Algorithm 
(NSGA)  and  Controlled  Elitist  Non-dominated  Sorting 
Genetic  Algorithms  (NSGA  II)  [13,14,15] 

After  a  new  population  is  formed,  the  strings  of  randomly 
selected  pairs  undergo  a  crossover  operation  with  a  probability 
Pc  and  mutation  with  a  probability  Pm.  We  have  found  that  Pc 
must  be  larger  than  0.7  and  Pm  should  be  smaller  than  0.1.  The 
size  of  the  population  should  be  larger  than  50.  Various 
crossover  types  can  be  selected  such  as  single  point,  two  point, 
and  arithmetic  crossover.  For  the  mutation  operation  also  we 
have  used  various  forms:  uniform  or  non-uniform  mutation. 
We  use  a  real  representation  for  the  gene  values.  A  detailed 
description  of  the  genetic  operators  is  given  in  reference  [8], 


2.4  Multiobjective  Optimization  with  Deterministic 
Algorithms 

These  objectives  allow  us  to  use  deterministic  gradient 
based  algorithms.  We  use  a  weighted  sum  approach  for  the 
multiobjective  optimization,  where  for  a  set  of  weights  for 
the  volume  and  surface  variance  we  perform  a  single 
objective  optimization  of 

fw^  fw  =wsfs  +  wvfv  (5) 

where  w$  ,  wv  >  0  are  the  surface  and  volume  importance 
factors,  respectively  and  ws  +  wy  =  1.  We  used  21 


2.6  Selecting  the  Solution  from  the  Pareto  Set 

After  the  last  generation  is  processed  by  the  SPEA,  FFGA, 
NPGA  or  NSGA  II  algorithm,  members  of  the  population  are 
expected  to  be  close  to  the  Pareto  frontier.  A  member  of  the 
non  dominated  set  is  selected  which  has  a  minimum  Euclidean 
distance  to  the  ideal  optimum.  The  ideal  point  is  defined  by 

the  minimum  values  (f"“n,  f"""  )  of  each  objective  function. 
The  distance  is  calculated  by  normalizing  each  objective  to  a 
maximum  value  of  1  using  the  corresponding  largest  objective 
value  found  in  the  population.  This  member  is  presented  as  the 


COIN 


solution  of  the  optimization  process.  Additionally  members 
are  selected  each  with  the  best  result  in  each  objective.  A  list 
is  produced  with  the  objective  values  for  all  the  members  of 
the  Pareto  set.  Additionally  the  user  can  examine  the  dose 
distributions  and  the  dose-volume  histogram  and  isodose 
contours  of  every  member  of  the  population.  Based  on  this 
information  of  the  trade-off  surface  of  the  various  objectives 

a  decision  maker  can 
select  the  best  result. 


Fig-2.  Flow  diagram  for  the  DVH- 
based  multiobjective  genetic  algorithm. 


Fig.3.  Comparison  of  the  COIN 
distributions  for  a  breast  implant  from  the 
multiobjective  genetic  algorithm  and  four 
conventional  single  objective  algorithms. 


3.  Results 

The  dose  variances 
are  calculated  from 
1000  -7-  4000  quasi- 
randomly  distributed 
sampling  points.  For 
the  COIN  based 
optimization  ~ 

100000  points  are 
generated.  The 
distances  of  these 
points  to  each  source 
dwell  position  r, 
more  precisely  the 
inverse  square 

distances  1/r2,  are 
stored  for  speed 
maximization  in 
look-up  tables.  We 
assume  a  invariant 
kernel  K(r)=  1/r  and 
ignore  any  spatial 
anisotropy,  namely 
attenuation  and 
scattering  effect. 
This  dosimetric 
simplification  has 
no  measurable 
influence  on  the 
results  of  the 
optimization. 

All  calculations 
presented  in  our 
study  have  been 
made  by  using  for 
the  mutation 

probability  Pm  a 
value  of  0.0065  and 
for  the  crossover 
probability  Pc  a 
value  of  0.85. 


The  flowchart  for  the  COIN  based  optimization  algorithm 
is  shown  in  Fig.  2.  For  each  member  of  the  population  for  a 
given  generation  a  renormalization  is  carried  out  according 
to  the  resulting  COIN  distribution,  so  that  the  maximum 
COIN  value  is  observed  at  D  =  Dref  [7].  The  dose 
prescription  is  realized  at  the  Dref,  the  isodose  value  resulting 
in  the  maximal  conformity.  This  results  generally  in  mean 
normalized  dose  values  at  the  surface  of  PTV  different  from 
1.0. 

The  multiobjective  genetic  algorithm,  which  uses  dose- 
volume  based  constraints,  produces  equivalent  or  even  better 
results  than  algorithms  which  were  based  on 


Fig.4.  (a)  Contours  of  a  rib  implant 
with  the  catheters  and  the  source 
dwell  positions.  On  the  PTV  surface 
sampling  points  are  shown  at  which 
the  dose  is  calculated  (b)  the  dose 
isosurface  obtained  from  the  dose 
optimization. 


phenomenological  methods  and  used  in  the  majority  of 
treatment  planning  systems  [16,17,18], 

As  an  example  in 
Fig.  3  the 

multiobjective  genetic 
algorithm  provides  a 
solution  with  a  more 
homogeneous  dose 
distribution  inside  the 
PTV  than  by 
conventional 
optimization 
algorithms  of  a 
treatment  planning 
system.  Due  to  the 
large  computational 
time  for  the  COIN 
based  optimization  we 
used  only  the  NPGA 
algorithm. 

An  example  of  the 
geometry  of  a  PTV  is 
shown  in  Fig.  5(a)Mic  including  the  catheters,  the  source 
dwell  positions  and  the  sampling  points  on  the  PTV  surface 
which  define  the  surface  variance.  In  Fig.  5(b)Mic  the 
isosurface  for  the  prescription  dose  is  shown,  which  should 
have  the  same  shape  as  the  PTV. 

The  deterministic  gradient  based  algorithms  are  very 
effective  in  generating  the  Pareto  front  using  a  summed 
weights  approach.  Powells  algorithm  which  does  not  require 
derivatives  is 
efficient  only 
for  implants 
with  a  small 
number  of 
sources  (time 
consuming), 
whereas  the 
gradient  based 
algorithms 
require  only 
1-2  minutes. 

Gradient 
based 

algorithms  are 
limited  by  the 
fact  that  they 
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Fig.5.  Pareto  front  obtained  by  the  gradient  based 
algorithm  and  with  the  SPEA  algorithm  with  and 
without  initialization. 


can  be  trapped  in  local  minima,  or  that  non  convex  regions  are 
not  accessible  using  the  weighted  sum  method  [19]. 

From  the  evolutionary  algorithms  NSGAII  and  SPEA  have 
been  found  to  produce  the  best  results,  since  it  applies  an 
elitism  and  sharing  mechanism.  For  implants  with  many 
sources  the  genetic  algorithms  used  converge  in  some  cases  to 
a  Pareto  set  which  was  far  away  from  the  true  Pareto  set.  Such 
an  example  for  an  implant  with  215  source  dwell  positions  is 
shown  in  Fig.  5.  The  SPEA  algorithm  converges  after  200 
generations  to  a  Pareto  front  which  is  very  small  and  far  from 
the  Pareto  set  generated  by  the  gradient  based  algorithms.  The 
optimization  path  is  shown  for  a  set  of  importance  factors  fv, 
fs  for  the  Polak-Ribiere  algorithm.  After  10  iterations  a  point 
on  the  Pareto  front  is  reached. 

Using  random  sets  of  decision  variables  we  have  found  for 
this  example  that  the  number  of  function  evaluations  required 
by  a  random  search  method  to  obtain  points  on  the  Pareto  front 


is  larger  than  1030  [8].  A  random  search  would  require  10 10 
times  more  function  evaluations  to  generate  points  on  the 
Pareto  set  found  by  the  SPEA  algorithm  without 
initialization.  Even  with  this  performance  the  SPEA 
algorithm  is  not  able  to  produce  points  on  the  Pareto  front 
found  by  the  deterministic  methods.  Using  a  few  members 
initialized  by  the  gradient  based  algorithm  the  multiobjective 
evolutionary  algorithms  reproduced  the  Pareto  fronts 
obtained  by  the  deterministic  algorithms.  Fig.  5.  For  a  more 
detailed  comparison  of  the  deterministic  and  evolutionary 
algorithms  see  reference  [8]. 

4.  Conclusions 

We  used  for  the  first  time  multiobjective  evolutionary 
anatomy  based  dose  optimization  algorithms  in  HDR 
brachytherapy  [6].  For  the  COfN-based  objectives  we  have 
found  that  multiobjective  evolutionary  algorithms  produced 
solutions  which  are  better  than  by  conventional  algorithms  in 
treatment  planning  systems  which  use  deterministic 
algorithms  and  catheter-oriented  objectives.  They  also  have 
the  problem  with  infeasible  negative  weights  which  they 
avoid  by  a  repair  mechanism  or  by  using  special  constraints 
to  the  objective  functions  in  order  to  reduce  their  numbers 
and  the  degree  of  the  violation. 

The  results  of  various  algorithms  for  the  variance  based 
objectives  have  been  compared  using  a  representative  set  of 
22  implants  encountered  in  clinical  practice.  We  have 
limited  our  study  to  cases  where  no  critical  structures  are 
considered.  Trade-off  surfaces  which  reveal  the  nature  of  the 
multiobjective  problem  of  the  dose  optimization  in 
brachytherapy  have  been  obtained.  Due  to  the  variety  of  the 
trade-off  surfaces  found,  which  depends  on  the  implant  and 
complex  catheter  geometry,  no  common  set  of  optimal 
importance  factors  exists.  Therefore  it  is  useful  to  determine 
the  Pareto  front  and  then  to  select  a  solution  according  to  its 
characteristics.  Pareto  sets  have  been  obtained  by  a 
deterministic  unconstrained  optimization  method  using  a 
simple  mapping  technique  which  transforms  the  linear  into  a 
quadratic  optimization  problem  and  removes  infeasible 
solutions  with  negative  dwell  position  weights.  The  gradient 
based  algorithms,  if  they  can  be  used,  are  very  effective 
because  they  converge  very  fast  and  generate  the  Pareto 
fronts  which  in  most  cases  are  much  better  than  the  Pareto 
front  obtained  by  evolutionary  multiobjective  algorithms. 

If  the  number  of  objectives  increases  then  the  number  of 
combinations  using  a  weighted  sum  approach  with 
deterministic  algorithms  increases.  Deterministic  methods 
are  not  efficient  for  non  analytic  complex  objectives  such  as 
used  by  the  COIN  based  method.  When  more  objectives  are 
included  then  a  non  convex  feasible  space  could  be  the  result 
[20],  A  combination  of  deterministic  and  evolutionary 
multiobjective  algorithms  seems  to  be  the  best  choice  for  a 
robust  and  efficient  multiobjective  dose  optimization  in 
HDR  brachytherapy. 

We  are  currently  studying  for  various  sets  of  objectives  the 
Pareto  fronts  using  multiobjective  evolutionary  algorithms 
and  if  possible  in  combination  with  deterministic  algorithms. 
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